
/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
/**
 * @author  John Miller
 * @version 1.0
 * @date    Sun Oct 25 18:41:26 EDT 2009
 * @see     LICENSE (MIT style license file).
 */

package scalation
package dynamics

import scala.math._
import advmath._
import util.Error

/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
/**
 * Given an unknown, time-dependent function y(t) governed by an Ordinary
 * Differential Equation (ODE) of the form y(t)' = f(t, y) where ' is d/dt,
 * compute y(t) using a 4th-order Runge-Kutta integrator.  Note: the integrateV
 * method for a system of ODE's is mixed in from the Integrator trait.
 */
object RungeKutta
       extends Integrator
{
    /*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
    /**
     * Compute y(t) governed by a differential equation using numerical integration
     * of the derivative function f(t, y) using a 4th-order Runge-Kutta method to
     * return the value of y(t) at time t.
     * @param f      the derivative function f(t, y)
     * @param y0     value of the y-function at time t0, y0 = y(t0)
     * @param t      the time value at which to compute y(t)
     * @param t0     the initial time
     * @param _step  the step size
     */
    def integrate (f: Derivative, y0: Double, t: Double,
                   t0: Double = 0., _step: Double = defaultStepSize): Double =
    {
        val t_t0       = t - t0                                     // time interval
   	val steps: Int = (round (t_t0 / _step)).asInstanceOf [Int]  // number of steps
   	var step       = t_t0 / steps.asInstanceOf [Double]         // adjusted step size
        var ti         = t0                                         // initialize ith time ti to t0
   	var y          = y0                                         // initialize y = f(t) to y0
        var k1 = 0.; var k2 = 0.; var k3 = 0.; var k4 = 0.

   	for (i <- 1 to steps) {
            ti += step                                       // take the next step
            if (ti > t) { step -= ti - t; ti = t }           // don't go past t

            k1  = step * f (ti, y)
            k2  = step * f (ti + step/2., y + k1/2.)
            k3  = step * f (ti + step/2., y + k2/2.)
            k4  = step * f (ti + step, y + k3)

            y  += (k1 + 2*k2 + 2*k3 + k4) / 6.

            if (abs (y) > Double.MaxValue / 10.) {
                flaw ("integrate", "probable overflow since y = " + y)
            } else if (abs (y) < Double.MinPositiveValue * 10.) {
                flaw ("integrate", "probable underflow since y = " + y)
            } // if
            if (i % 1000 == 0) println ("integrate: iteration " + i + " ti = " + ti + " y = " + y)
   	} // for

        y         // the value of the function at time t, y = f(t)
    } // integrate

} // RungeKutta object

/*:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::*/
/**
 * This object is used to test the RungeKutta object.
 */
object RungeKuttaTest extends App
{
    import scalation.dynamics.RungeKutta._

    val y0 = 1.
    val t  = 2.

    def derv1 (t: Double, y: Double) = 2. * t  // solution to differential equation is t^2
    println ("\n==> at t = " + t + " y = " + integrate (derv1, y0, t))
    println ("\n==> t^2 + c = " + 5)

    def derv2 (t: Double, y: Double) = y       // solution to differential equation is e^t
    println ("\n==> at t = " + t + " y = " + integrate (derv2, y0, t))
    println ("\n==> e = " + E + " e^2 = " + E * E)

    def derv3 (t: Double, y: Double) = t + y   // solution to differential equation is ?
    println ("\n==> at t = " + t + " y = " + integrate (derv3, y0, t))

    println ("\n==> at t = " + t + " y = " + 
             integrateV (Array (derv1, derv2), Vec(1., 2.), t))

} // RungeKuttaTest object

